/*
 * NumUtils.cpp
 *
 *  Created on: 2014-5-28
 *      Author: cpp
 */

#include "NumUtils.h"
#include <strstream>
#include<math.h>
using namespace std;
double NumUtils::MAXLOG = 7.09782712893383996732E2;
double NumUtils::SQRTH = 7.07106781186547524401E-1;
string NumUtils::ftos(double value) {
	strstream ss;
	ss << value;
	string result = "";
	ss >> result;
	return result;
}

int NumUtils::sgn(double x) {
	if (x > 0) {
		return 1;
	} else if(x==0){
		return 0;
	}else{
		return -1;
	}
}

double NumUtils::copysign(double x, double y) {
	return x * sgn(y);
}

const double NumUtils::U = 0.5 * pow(FLT_RADIX, (1 - DBL_MANT_DIG));
const double NumUtils::e = 1.0e-16;
const double NumUtils::ERR = 0.0001;

/**
 * Returns the area under the Normal (Gaussian) probability density function,
 * integrated from minus infinity to <tt>x</tt> (assumes mean is zero,
 * variance is one).
 *
 * <pre>
 *                            x
 *                             -
 *                   1        | |          2
 *  normal(x)  = ---------    |    exp( - t /2 ) dt
 *               sqrt(2pi)  | |
 *                           -
 *                          -inf.
 *
 *             =  ( 1 + erf(z) ) / 2
 *             =  erfc(z) / 2
 * </pre>
 *
 * where <tt>z = x/sqrt(2)</tt>. Computation is via the functions
 * <tt>errorFunction</tt> and <tt>errorFunctionComplement</tt>.
 *
 * @param a the z-value
 * @return the probability of the z value according to the normal pdf
 */
double NumUtils::normalProbability(double a) {

	double x, y, z;

	x = a * SQRTH;
	z = fabs(x);

	if (z < SQRTH) {
		//y = 0.5 + 0.5 * errorFunction(x);
		y = 0.5 + 0.5 * erf(x);
	} else {
		//y = 0.5 * errorFunctionComplemented(z);
		y = 0.5 *erfc(x);;//;/
		if (x > 0) {
			y = 1.0 - y;
		}
	}
	return y;
}

double NumUtils::digamma(double x){
    double p;  
    x=x+6;  
    p=1/(x*x);  
    p=(((0.004166666666667*p-0.003968253986254)*p+  
    0.008333333333333)*p-0.083333333333333)*p;  
    p=p+log(x)-0.5/x-1/(x-1)-1/(x-2)-1/(x-3)-1/(x-4)-1/(x-5)-1/(x-6);  
    return p;
}

double NumUtils::trigamma(double x){
    double p;  
    int i;   
    x=x+6;  
    p=1/(x*x);  
    p=(((((0.075757575757576*p-0.033333333333333)*p+0.0238095238095238)  
         *p-0.033333333333333)*p+0.166666666666667)*p+1)/x+0.5*p;  
    for (i=0; i<6 ;i++)  
    {  
        x=x-1;  
        p=1/(x*x)+p;  
    }
    return(p);
}
